1 Packages
library(ggplot2) #for plotting
library(tidyverse)## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ tibble 3.0.1 ✓ dplyr 1.0.0
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ✓ purrr 0.3.4## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(dplyr)
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine2 Read in Data
The Pinery dataset was compiled in excel from monthly mean temperature and total precipitation data using open access Environment Canada records for stations within 30km of the Pinery with at least 10 consecutive years of data (http://climate.weather.gc.ca/historical_data/search_historic_data_e.html).
The London datasets are homogenized data from Environment Canada that were converted to tidy format in excel.
pinery <- read.csv("data/pinery_climate.csv")
pinery$station <- as.factor(pinery$station) #Make a factor
str(pinery)## 'data.frame': 2931 obs. of 5 variables:
## $ year : int 1882 1882 1882 1882 1882 1882 1883 1883 1883 1883 ...
## $ month : int 7 8 9 10 11 12 4 5 6 7 ...
## $ mean_temp : num NA NA 16.1 11.4 2.9 -3.9 NA NA NA NA ...
## $ total_precip: num 15.5 182.6 57.9 41.7 28.7 ...
## $ station : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 8 8 8 8 ...london_precip <- read.csv("data/london_precip.csv")
london_temp <- read.csv("data/london_temp.csv")
london_precip$month <- as.factor(london_precip$month) #Make a factor
london_temp$month <- as.factor(london_temp$month) #Make a factor climateNA <- read_csv("data/climateNA_monthly.csv") ## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_double(),
## Year = col_double(),
## Month = col_double(),
## Tmax = col_double(),
## Tmin = col_double(),
## Tave = col_double(),
## PPT = col_double(),
## Rad = col_double(),
## Eref = col_double(),
## CMD = col_double(),
## RH = col_double()
## )climateNA.temp <- climateNA %>% select(Year, Month, Tave)
climateNA.precip <- climateNA %>% select(Year, Month, PPT)3 Create Regional Data
3.1 Clean Data
Here I go through each station for the temp and precip data, removing years where all 12 months aren’t present so that I can later calculate the mean values for each year. I did this manually a while ago and there is probably a more efficient way to do it, however it works and it’s already done so I left it as is.
#Arkona
Arkona <- pinery %>% filter(station == "Arkona") #Filter out station of interest
Arkona_T <- Arkona %>% select(-total_precip) %>% filter(!is.na(Arkona$mean_temp)==T) #remove precipitation data and "NA"
table(Arkona_T$year) #Generate table after cleaning##
## 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897
## 4 11 11 12 12 12 12 12 12 12 12 5 12 12 12 12
## 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913
## 12 12 10 12 4 12 12 12 12 12 12 12 12 12 12 12
## 1914 1915
## 12 3Arkona_T2 <- Arkona_T %>% filter(year != 1882 & year != 1883 &
year != 1884 & year != 1893 & year != 1900 & year != 1902 & year != 1915) #remove years with < 12 months
table(Arkona_T2$year) #Generate table after cleaning##
## 1885 1886 1887 1888 1889 1890 1891 1892 1894 1895 1896 1897 1898 1899 1901 1903
## 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
## 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914
## 12 12 12 12 12 12 12 12 12 12 12Arkona_P <- Arkona %>% select(-mean_temp) %>% filter(!is.na(Arkona$total_precip)==T) #remove temperature data and "NA"
#table(Arkona_P$year) #show how many observations in each year
Arkona_P2 <- Arkona_P %>% filter(year != 1882 & year != 1893 & year != 1902 & year != 1915)
#table(Arkona_P2$year) #Generate table after cleaning
#Centralia A
CentraliaA <- pinery %>% filter(station == "Centralia A")
CentraliaA_T <- CentraliaA %>% select(-total_precip) %>% filter(!is.na(CentraliaA$mean_temp)==T) #remove "NA"
#table(CentraliaA_T$year)
CentraliaA_T2 <- CentraliaA_T %>% filter(year != 1945 & year != 1946 & year != 1947 & year != 1966)
#table(CentraliaA_T2$year)
CentraliaA_P <- CentraliaA %>% select(-mean_temp) %>% filter(!is.na(CentraliaA$total_precip)==T) #remove "NA"
#table(CentraliaA_P$year)
CentraliaA_P2 <- CentraliaA_P %>% filter(year != 1946 & year != 1947 & year != 1966)
#table(CentraliaA_P2$year)
#Exeter
Exeter <- pinery %>% filter(station == "Exeter")
Exeter_T <- Exeter %>% select(-total_precip) %>% filter(!is.na(Exeter$mean_temp)==T) #remove "NA"
#table(Exeter_T$year)
Exeter_T2 <- Exeter_T %>% filter(year != 1967 & year != 1976 & year != 1985)
#table(Exeter_T2$year)
Exeter_P <- Exeter %>% select(-mean_temp) %>% filter(!is.na(Exeter$total_precip)==T) #remove "NA"
#table(Exeter_P$year)
Exeter_P2 <- Exeter_P %>% filter(year != 1961 & year != 1962 & year != 1964 & year != 1976 & year != 1985)
#table(Exeter_P2$year)
#Forest
Forest <- pinery %>% filter(station == "Forest")
Forest_T <- Forest %>% select(-total_precip) %>% filter(!is.na(Forest$mean_temp)==T) #remove "NA"
#table(Forest_T$year)
Forest_T2 <- Forest_T %>% filter(year != 1924 & year != 1940 & year != 1962 & year != 1963)
#table(Forest_T2$year)
Forest_P <- Forest %>% select(-mean_temp) %>% filter(!is.na(Forest$total_precip)==T) #remove "NA"
#table(Forest_P$year)
Forest_P2 <- Forest_P %>% filter(year > 1933 & year != 1937 & year!= 1943 & year != 1955 & year != 1960 & year < 1960)
#table(Forest_P2$year)
#Lucan
Lucan <- pinery %>% filter(station == "Lucan")
Lucan_T <- Lucan %>% select(-total_precip) %>% filter(!is.na(Lucan$mean_temp)==T) #remove "NA"
#table(Lucan_T$year)
Lucan_T2 <- Lucan_T %>% filter(year != 1915 & year != 1924 & year != 1948 & year != 1949 & year != 1950
& year != 1951 & year != 1952 & year != 1953 & year != 1954 & year != 1955
& year != 1956 & year != 1957 & year != 1958 & year != 1959 & year != 1961)
#table(Lucan_T2$year)
Lucan_P <- Lucan %>% select(-mean_temp) %>% filter(!is.na(Lucan$total_precip)==T) #remove "NA"
#table(Lucan_P$year)
Lucan_P2 <- Lucan_P %>% filter(year != 1915 & year != 1923 & year != 1943 & year < 1948)
#table(Lucan_P2$year)
#Pinery
Pinery <- pinery %>% filter(station == "Pinery")
Pinery_T <- Pinery %>% select(-total_precip) %>% filter(!is.na(Pinery$mean_temp)==T) #remove "NA"
#table(Pinery_T$year)
Pinery_T2 <- Pinery_T %>% filter(year != 1979 & year != 1984)
#table(Pinery_T2$year)
Pinery_P <- Pinery %>% select(-mean_temp) %>% filter(!is.na(Pinery$total_precip)==T) #remove "NA"
#table(Pinery_P$year)
Pinery_P2 <- Pinery_P %>% filter(year != 1979 & year != 1984)
#table(Pinery_P2$year)
#Strathrow-Mullifary
Strathroy <- pinery %>% filter(station == "Strathroy-Mullifary")
Strathroy_T <- Strathroy %>% select(-total_precip) %>% filter(!is.na(Strathroy$mean_temp)==T) #remove "NA"
#table(Strathroy_T$year)
Strathroy_T2 <- Strathroy_T %>% filter(year != 2009 & year != 2010 & year != 2013 & year != 2014 & year != 2015 & year < 2017)
#table(Strathroy_T2$year)
Strathroy_P <- Strathroy %>% select(-mean_temp) %>% filter(!is.na(Strathroy$total_precip)==T) #remove "NA"
#table(Strathroy_P$year)
Strathroy_P2 <- Strathroy_P %>% filter(year != 2009 & year != 2010 & year != 2012 & year != 2013
& year != 2014 & year != 2015 & year < 2017)
#table(Strathroy_P2$year)
#Thedford (no temp)
Thedford <- pinery %>% filter(station == "Thedford")
Thedford_P <- Thedford %>% select(-mean_temp) %>% filter(!is.na(Thedford$total_precip)==T) #remove "NA"
#table(Thedford_P$year)
Thedford_P2 <- Thedford_P %>% filter(year != 1883 & year != 1897)
#table(Thedford_P2$year)
#Watford (can't use precipitation data)
Watford <- pinery %>% filter(station == "Watford")
Watford_T <- Watford %>% select(-total_precip) %>% filter(!is.na(Watford$mean_temp)==T) #remove "NA"
#table(Watford_T$year)
Watford_T2 <- Watford_T %>% filter(year == 1927 | year == 1928)
#table(Watford_T2$year)
Watford_P <- Watford %>% select(-mean_temp) %>% filter(!is.na(Watford$total_precip)==T) #remove "NA"
#table(Watford_P$year) Re-combine Data
pinery_temp <- rbind(Arkona_T2, CentraliaA_T2, Exeter_T2, Forest_T2, Lucan_T2, Pinery_T2, Strathroy_T2, Watford_T2)
pinery_precip <- rbind(Arkona_P2, CentraliaA_P2, Exeter_P2, Forest_P2, Lucan_P2, Pinery_P2, Strathroy_P2, Thedford_P2)Ensure the cleaning worked by plotting histograms of the number of observations (months) in each year (should only be multiples of 12)
#histogram (temperature)
graph <- ggplot(pinery_temp, aes(x = year)) +
geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
geom_hline(yintercept = 24, linetype="solid", color = "red") +
geom_hline(yintercept = 36, linetype="solid", color = "red") +
theme_bw()
graph
#histogram (precipitation)
graph <- ggplot(pinery_precip, aes(x = year)) +
geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
geom_hline(yintercept = 24, linetype="solid", color = "red") +
theme_bw()
graph
3.2 Create Composite Record
Calculate one mean temperature and total precipitation for each year/month where more than one station of data is available
#Temperature
str(pinery_temp)## 'data.frame': 1944 obs. of 4 variables:
## $ year : int 1885 1885 1885 1885 1885 1885 1885 1885 1885 1885 ...
## $ month : int 1 2 3 4 5 6 7 8 9 10 ...
## $ mean_temp: num -9.4 -13.5 -8.9 4.2 12.1 16.4 21.1 17 14.8 8 ...
## $ station : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 1 1 1 1 ...pinery_temp_mean <- pinery_temp %>%
group_by(year,month) %>%
summarise(mean_temp_year = mean(mean_temp), sd_temp = sd(mean_temp), n = n()) %>%
mutate(se_temp = sd_temp / sqrt(n), lower = (mean_temp_year - se_temp), upper = (mean_temp_year + se_temp)) %>% round(2)## `summarise()` regrouping output by 'year' (override with `.groups` argument)pinery_temp_mean$month <- as.factor(pinery_temp_mean$month)
#Precipitation
pinery_precip_mean <- pinery_precip %>%
group_by(year,month) %>%
summarise(total_precip_year = mean(total_precip), sd_precip = sd(total_precip), n = n()) %>%
mutate(se_precip = sd_precip / sqrt(n), lower = total_precip_year - sd_precip, upper = total_precip_year + sd_precip) %>% round(2)## `summarise()` regrouping output by 'year' (override with `.groups` argument)head(pinery_precip_mean)## # A tibble: 6 x 8
## # Groups: year [1]
## year month total_precip_year sd_precip n se_precip lower upper
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1883 1 30.2 NA 1 NA NA NA
## 2 1883 2 94.7 NA 1 NA NA NA
## 3 1883 3 70.6 NA 1 NA NA NA
## 4 1883 4 33 NA 1 NA NA NA
## 5 1883 5 112. NA 1 NA NA NA
## 6 1883 6 158. NA 1 NA NA NAggplot(pinery_precip_mean, aes(x = year)) +
geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey")
ggplot(pinery_temp_mean, aes(x = year)) +
geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey")
3.3 Fill Empty Years
3.3.1 Temperature
Fill empty months-years with data from London
p.t <- pinery_temp_mean %>% select(year, month, mean_temp_year)
l.t <- london_temp %>% filter(month == "Jan" |
month == "Feb" |
month == "Mar" |
month == "Apr" |
month == "May" |
month == "Jun" |
month == "Jul" |
month == "Aug" |
month == "Sep" |
month == "Oct" |
month == "Nov" |
month == "Dec") %>%
mutate(month = recode(month, "Jan"="1",
"Feb"="2",
"Mar"="3",
"Apr"="4",
"May"="5",
"Jun"="6",
"Jul"="7",
"Aug"="8",
"Sep"="9",
"Oct"="10",
"Nov"="11",
"Dec"="12"))PL.temp <- full_join(l.t, p.t) %>%
rename("pinery" = "mean_temp_year") %>%
rename("london" = "mean_temp")## Joining, by = c("year", "month")head(PL.temp)## year month london pinery
## 1 1883 1 NA NA
## 2 1884 1 -10.1 NA
## 3 1885 1 -9.2 -9.4
## 4 1886 1 -7.2 -7.2
## 5 1887 1 -8.3 -8.4
## 6 1888 1 -10.1 -9.6Select all rows with regional data - source = Pinery Select all rows with NA - source = London
PL.1 <- PL.temp %>% #use london value where pinery is absent
filter(is.na(pinery)) %>%
mutate(temp = london)
PL.1$source <- "London"
PL.2 <- PL.temp %>% # use pinery value where pinery is present
filter(!is.na(pinery)) %>%
mutate(temp = pinery)
PL.2$source <- "Pinery"
PL.temp2 <- rbind(PL.1, PL.2) %>%
arrange(year)
(nrow(PL.2) + nrow(PL.1)) == nrow(PL.temp) #check equal observations## [1] TRUEPL.temp2 %>% filter(is.na(temp))## year month london pinery temp source
## 1 1883 1 NA NA NA London
## 2 1883 2 NA NA NA London
## 3 2014 9 NA NA NA LondonFill leftover NAs with climateNA values
head(climateNA.temp)## # A tibble: 6 x 3
## Year Month Tave
## <dbl> <dbl> <dbl>
## 1 1901 1 -4.2
## 2 1902 1 -5.6
## 3 1903 1 -5.1
## 4 1904 1 -9.5
## 5 1905 1 -7.5
## 6 1906 1 -0.7climateNA.temp <- climateNA.temp %>%
rename("month" = "Month") %>%
rename("year" = "Year") %>%
rename("climateNA" = "Tave")
climateNA.temp$month <- as.factor(climateNA.temp$month)
PL.temp3 <- full_join(PL.temp2, climateNA.temp)## Joining, by = c("year", "month")PL.temp3 <- PL.temp3[, c(1:4,7,5,6)] #re-order columns
head(PL.temp3)## year month london pinery climateNA temp source
## 1 1883 1 NA NA NA NA London
## 2 1883 2 NA NA NA NA London
## 3 1883 3 -6.5 NA NA -6.5 London
## 4 1883 4 5.1 NA NA 5.1 London
## 5 1883 5 10.0 NA NA 10.0 London
## 6 1883 6 17.3 NA NA 17.3 LondonPL.3 <- PL.temp3 %>%
filter(is.na(temp)) %>%
mutate(temp = climateNA)
PL.3$source <- "ClimateNA"
PL.4 <- PL.temp3 %>%
filter(!is.na(temp))
PL.temp4 <- rbind(PL.3, PL.4) %>%
arrange(year) %>%
filter(year > 1890) #missing values from 1890s remain
(nrow(PL.3) + nrow(PL.4)) == nrow(PL.temp3) #check equal observations## [1] TRUEPL.temp4 %>% filter(is.na(temp)) #values still missing## [1] year month london pinery climateNA temp source
## <0 rows> (or 0-length row.names)ggplot(PL.temp4, aes(x = year)) +
geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
theme_bw()
3.3.2 Precipitation
p.p <- pinery_precip_mean %>% select(year, month, total_precip_year)
l.p <- london_precip %>% filter(month == "Jan" |
month == "Feb" |
month == "Mar" |
month == "Apr" |
month == "May" |
month == "Jun" |
month == "Jul" |
month == "Aug" |
month == "Sep" |
month == "Oct" |
month == "Nov" |
month == "Dec") %>%
mutate(month = recode(month, "Jan"="1",
"Feb"="2",
"Mar"="3",
"Apr"="4",
"May"="5",
"Jun"="6",
"Jul"="7",
"Aug"="8",
"Sep"="9",
"Oct"="10",
"Nov"="11",
"Dec"="12"))p.p$month <- as.factor(p.p$month)
PL.precip <- full_join(l.p, p.p) %>%
rename("pinery" = "total_precip_year") %>%
rename("london" = "total_precip")## Joining, by = c("year", "month")PL.5 <- PL.precip %>% #use london value where pinery is absent
filter(is.na(pinery)) %>%
mutate(precip = london)
PL.5$source <- "London"
PL.6 <- PL.precip %>% # use pinery value where pinery is present
filter(!is.na(pinery)) %>%
mutate(precip = pinery)
PL.6$source <- "Pinery"
PL.precip2 <- rbind(PL.5, PL.6) %>%
arrange(year)
(nrow(PL.5) + nrow(PL.6)) == nrow(PL.precip2) #check equal observations## [1] TRUEPL.precip2 %>% filter(is.na(precip)) #values still missing## year month london pinery precip source
## 1 2009 4 NA NA NA London
## 2 2009 5 NA NA NA London
## 3 2009 6 NA NA NA London
## 4 2009 7 NA NA NA London
## 5 2009 8 NA NA NA London
## 6 2009 9 NA NA NA London
## 7 2009 10 NA NA NA London
## 8 2010 4 NA NA NA London
## 9 2010 5 NA NA NA London
## 10 2010 6 NA NA NA London
## 11 2010 7 NA NA NA London
## 12 2010 8 NA NA NA London
## 13 2010 9 NA NA NA London
## 14 2010 10 NA NA NA London
## 15 2012 4 NA NA NA London
## 16 2012 5 NA NA NA London
## 17 2012 6 NA NA NA London
## 18 2012 7 NA NA NA London
## 19 2012 8 NA NA NA London
## 20 2012 9 NA NA NA London
## 21 2012 10 NA NA NA London
## 22 2013 4 NA NA NA London
## 23 2013 5 NA NA NA London
## 24 2013 6 NA NA NA London
## 25 2013 7 NA NA NA London
## 26 2013 8 NA NA NA London
## 27 2013 9 NA NA NA London
## 28 2013 10 NA NA NA London
## 29 2014 4 NA NA NA London
## 30 2014 5 NA NA NA London
## 31 2014 6 NA NA NA London
## 32 2014 7 NA NA NA London
## 33 2014 8 NA NA NA London
## 34 2014 9 NA NA NA London
## 35 2014 10 NA NA NA London
## 36 2015 4 NA NA NA London
## 37 2015 5 NA NA NA London
## 38 2015 6 NA NA NA London
## 39 2015 7 NA NA NA London
## 40 2015 8 NA NA NA London
## 41 2015 9 NA NA NA London
## 42 2015 10 NA NA NA London
## 43 2017 4 NA NA NA London
## 44 2017 5 NA NA NA London
## 45 2017 6 NA NA NA London
## 46 2017 7 NA NA NA London
## 47 2017 8 NA NA NA London
## 48 2017 9 NA NA NA London
## 49 2017 10 NA NA NA London
## 50 2017 11 NA NA NA London
## 51 2017 12 NA NA NA LondonFill leftover NAs with climateNA values
head(climateNA.temp)## # A tibble: 6 x 3
## year month climateNA
## <dbl> <fct> <dbl>
## 1 1901 1 -4.2
## 2 1902 1 -5.6
## 3 1903 1 -5.1
## 4 1904 1 -9.5
## 5 1905 1 -7.5
## 6 1906 1 -0.7climateNA.precip <- climateNA.precip %>%
rename("month" = "Month") %>%
rename("year" = "Year") %>%
rename("climateNA" = "PPT")
climateNA.precip$month <- as.factor(climateNA.precip$month)
PL.precip3 <- full_join(PL.precip2, climateNA.precip)## Joining, by = c("year", "month")PL.precip3 <- PL.precip3[, c(1:4,7,5,6)] #re-order columns
head(PL.temp3)## year month london pinery climateNA temp source
## 1 1883 1 NA NA NA NA London
## 2 1883 2 NA NA NA NA London
## 3 1883 3 -6.5 NA NA -6.5 London
## 4 1883 4 5.1 NA NA 5.1 London
## 5 1883 5 10.0 NA NA 10.0 London
## 6 1883 6 17.3 NA NA 17.3 LondonPL.7 <- PL.precip3 %>%
filter(is.na(precip)) %>%
mutate(precip = climateNA)
PL.7$source <- "ClimateNA"
PL.8 <- PL.precip3 %>%
filter(!is.na(precip))
PL.7## year month london pinery climateNA precip source
## 1 2009 4 NA NA 132 132 ClimateNA
## 2 2009 5 NA NA 67 67 ClimateNA
## 3 2009 6 NA NA 108 108 ClimateNA
## 4 2009 7 NA NA 72 72 ClimateNA
## 5 2009 8 NA NA 90 90 ClimateNA
## 6 2009 9 NA NA 40 40 ClimateNA
## 7 2009 10 NA NA 116 116 ClimateNA
## 8 2010 4 NA NA 64 64 ClimateNA
## 9 2010 5 NA NA 103 103 ClimateNA
## 10 2010 6 NA NA 132 132 ClimateNA
## 11 2010 7 NA NA 100 100 ClimateNA
## 12 2010 8 NA NA 36 36 ClimateNA
## 13 2010 9 NA NA 76 76 ClimateNA
## 14 2010 10 NA NA 57 57 ClimateNA
## 15 2012 4 NA NA 58 58 ClimateNA
## 16 2012 5 NA NA 50 50 ClimateNA
## 17 2012 6 NA NA 52 52 ClimateNA
## 18 2012 7 NA NA 67 67 ClimateNA
## 19 2012 8 NA NA 70 70 ClimateNA
## 20 2012 9 NA NA 73 73 ClimateNA
## 21 2012 10 NA NA 128 128 ClimateNA
## 22 2013 4 NA NA 115 115 ClimateNA
## 23 2013 5 NA NA 71 71 ClimateNA
## 24 2013 6 NA NA 131 131 ClimateNA
## 25 2013 7 NA NA 91 91 ClimateNA
## 26 2013 8 NA NA 85 85 ClimateNA
## 27 2013 9 NA NA 51 51 ClimateNA
## 28 2013 10 NA NA 123 123 ClimateNA
## 29 2014 4 NA NA 81 81 ClimateNA
## 30 2014 5 NA NA 81 81 ClimateNA
## 31 2014 6 NA NA 88 88 ClimateNA
## 32 2014 7 NA NA 90 90 ClimateNA
## 33 2014 8 NA NA 103 103 ClimateNA
## 34 2014 9 NA NA 95 95 ClimateNA
## 35 2014 10 NA NA 84 84 ClimateNA
## 36 2015 4 NA NA 74 74 ClimateNA
## 37 2015 5 NA NA 95 95 ClimateNA
## 38 2015 6 NA NA 118 118 ClimateNA
## 39 2015 7 NA NA 65 65 ClimateNA
## 40 2015 8 NA NA 79 79 ClimateNA
## 41 2015 9 NA NA 59 59 ClimateNA
## 42 2015 10 NA NA 68 68 ClimateNA
## 43 2017 4 NA NA 118 118 ClimateNA
## 44 2017 5 NA NA 122 122 ClimateNA
## 45 2017 6 NA NA 116 116 ClimateNA
## 46 2017 7 NA NA 75 75 ClimateNA
## 47 2017 8 NA NA 82 82 ClimateNA
## 48 2017 9 NA NA 86 86 ClimateNA
## 49 2017 10 NA NA 107 107 ClimateNA
## 50 2017 11 NA NA 110 110 ClimateNA
## 51 2017 12 NA NA 50 50 ClimateNA
## 52 2018 1 NA NA 66 66 ClimateNA
## 53 2018 2 NA NA 107 107 ClimateNA
## 54 2018 3 NA NA 68 68 ClimateNA
## 55 2018 4 NA NA 102 102 ClimateNA
## 56 2018 5 NA NA 67 67 ClimateNA
## 57 2018 6 NA NA 97 97 ClimateNA
## 58 2018 7 NA NA 89 89 ClimateNA
## 59 2018 8 NA NA 180 180 ClimateNA
## 60 2018 9 NA NA 71 71 ClimateNA
## 61 2018 10 NA NA 100 100 ClimateNA
## 62 2018 11 NA NA 136 136 ClimateNA
## 63 2018 12 NA NA 58 58 ClimateNAPL.precip4 <- rbind(PL.7, PL.8) %>%
arrange(year) %>%
filter(year > 1890)
(nrow(PL.7) + nrow(PL.8)) == nrow(PL.precip3) #check equal observations## [1] TRUEPL.precip4 %>% filter(is.na(precip)) #values still missing## [1] year month london pinery climateNA precip source
## <0 rows> (or 0-length row.names)ggplot(PL.precip4, aes(x = year)) +
geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
theme_bw()
3.3.3 Combine
Fix levels
PL.temp4 <- droplevels(PL.temp4) #drops un-used levels
PL.temp4$month <- factor(PL.temp4$month, levels = c("1","2","3","4", "5", "6", "7", "8", "9", "10", "11", "12")) #re-level
levels(PL.temp4$month)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12"PL.precip4 <- droplevels(PL.precip4) #drops un-used levels
PL.precip4$month <- factor(PL.precip4$month, levels = c("1","2","3","4", "5", "6", "7", "8", "9", "10", "11", "12")) #re-level
levels(PL.precip4$month)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12"PL.temp.only <- PL.temp4 %>% select(year, month, temp)
PL.precip.only <- PL.precip4 %>% select(year, month, precip)
PL.clean <- full_join(PL.temp.only, PL.precip.only, by = c("year", "month"))
str(PL.clean)## 'data.frame': 1536 obs. of 4 variables:
## $ year : num 1891 1891 1891 1891 1891 ...
## $ month : Factor w/ 12 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ temp : num -4.2 -1.9 -2.5 7 10.9 18.8 17.5 18.8 16.9 8.6 ...
## $ precip: num 53.6 72.5 89 45.4 40 ...PL.clean %>% filter(is.na(temp)) #values still missing## [1] year month temp precip
## <0 rows> (or 0-length row.names)nrow(PL.clean) == (2019-1891)*12 #check 12 months / year for the 128 years## [1] TRUE4 Monthly Data
4.1 Climograph
study_t <- PL.clean %>%
group_by(month) %>%
summarize(mean_temp = mean(temp) %>% round(2)) ## `summarise()` ungrouping output (override with `.groups` argument)head(study_t)## # A tibble: 6 x 2
## month mean_temp
## <fct> <dbl>
## 1 1 -5.68
## 2 2 -5.68
## 3 3 -0.46
## 4 4 6.38
## 5 5 12.8
## 6 6 18.2study_p <- PL.clean %>%
group_by(month) %>%
summarize(total_precip= mean(precip), sd_precip = sd(precip), n = n()) %>%
mutate(se_precip = sd_precip / sqrt(n), lower = total_precip - se_precip, upper = total_precip + se_precip %>% round(2)) ## `summarise()` ungrouping output (override with `.groups` argument)head(study_p)## # A tibble: 6 x 7
## month total_precip sd_precip n se_precip lower upper
## <fct> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 1 82.8 33.3 128 2.94 79.8 85.7
## 2 2 67.3 28.5 128 2.52 64.8 69.8
## 3 3 65.7 28.3 128 2.50 63.2 68.2
## 4 4 74.6 32.0 128 2.83 71.8 77.5
## 5 5 81.5 36.9 128 3.26 78.2 84.8
## 6 6 81.0 37.7 128 3.33 77.6 84.3p.climograph <- ggplot() +
geom_boxplot(data = PL.clean, aes(x = month, y = temp), fill = "red3", alpha = 0.7) +
geom_path(data = study_p, aes(x = month, y = total_precip-75, group=1), col= "steelblue4", size = 1) +
geom_ribbon(data = study_p, aes(x = month, ymin = (lower-75), ymax = (upper-75), group=1), fill = "steelblue4", alpha = 0.2) +
scale_y_continuous("Temperature (°C)", sec.axis = sec_axis(~ . + 75, name = "Precipitation (mm)")) +
scale_x_discrete("Month") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p.climograph
#ggsave('p.climograph.pdf', p.climograph, units = 'cm', width = 20, height = 12)ggplot(data = PL.temp4, aes(x=year, y = temp)) +
geom_point(aes(col=source), size=0.3) +
geom_line(alpha = 0.6) +
facet_grid(month~.) +
ylab("Mean Monthly Temperature (°C)") +
theme_classic() 
ggplot(data = PL.temp4, aes(x=year, y = temp)) +
geom_point(aes(col=source), size=0.3) +
facet_grid(month~source) +
ylab("Mean Monthly Temperature (°C)") +
theme_classic() 
ggplot(data = PL.precip4, aes(x=year, y = precip)) +
geom_point(aes(col=source), size=0.3) +
geom_line(alpha = 0.6) +
facet_grid(month~.) +
ylab("Total Monthly Precipitation (mm)") +
theme_classic() 
ggplot(data = PL.precip4, aes(x=year, y = precip)) +
geom_point(aes(col=source), size=0.3) +
facet_grid(month~source) +
ylab("Total Monthly Precipitation (mm)") +
theme_classic() 
ggplot(data = PL.temp4, aes(x = month, y = temp)) +
geom_boxplot(alpha = 0.7, outlier.shape=NA) +
geom_jitter(aes(col=source), size=0.6, alpha=0.5) +
ylab("Mean Monthly Temperature (°C)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
ggplot(data = PL.precip4, aes(x = month, y = precip)) +
geom_boxplot(alpha = 0.7, outlier.shape=NA) +
geom_jitter(aes(col=source), size=0.6, alpha=0.5) +
ylab("Total Monthly Precipitation (mm)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
5 Annual Data
5.1 Clean Record
clean_temp_annual <- PL.clean %>% group_by(year) %>% summarise(mean_temp = mean(temp))## `summarise()` ungrouping output (override with `.groups` argument)clean_precip_annual <- PL.clean %>% group_by(year) %>% summarise(mean_precip = mean(precip))## `summarise()` ungrouping output (override with `.groups` argument)ggplot() +
geom_line(data=clean_temp_annual, aes(x = year, y = mean_temp), colour="lightsalmon4") +
#Aesthetics
xlab("Year") +
ylab("Mean Annual Temperature (°C)") +
scale_x_continuous(limits = c(1890, 2020), breaks=seq(1900,2020,20)) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
ggplot() +
geom_line(data=clean_precip_annual, aes(x = year, y = mean_precip), colour="lightsalmon4") +
xlab("Year") +
ylab("Total Annual Precipitation (mm)") +
scale_x_continuous(limits = c(1890, 2020), breaks=seq(1900,2020,20)) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
5.2 Individual Records
london_temp_annual <- london_temp %>% filter(month == "Annual")
london_precip_annual <- london_precip %>% filter(month == "Annual")
pinery_temp_annual <- pinery_temp_mean %>% group_by(year) %>% summarise(mean_temp = mean(mean_temp_year))## `summarise()` ungrouping output (override with `.groups` argument)pinery_precip_annual <- pinery_precip_mean %>% group_by(year) %>% summarise(total_precip = sum(total_precip_year))## `summarise()` ungrouping output (override with `.groups` argument)london_temp_2 <- london_temp_annual %>% select(-month)
london_temp_2$location <- "London"
pinery_temp_annual$location <- "Pinery"
annual_temp <- rbind(london_temp_2,pinery_temp_annual) #join data sets
london_precip_2 <- london_precip_annual %>% select(-month)
london_precip_2$location <- "London"
pinery_precip_annual$location <- "Pinery"
annual_precip <- rbind(london_precip_2,pinery_precip_annual) #join data setsannual_temp_plot <- ggplot() +
# London Data
geom_line(data=london_temp_annual, aes(x = year, y = mean_temp), colour = "lightseagreen", linetype = "dashed", alpha = 0.8) +
# Pinery Data
geom_line(data=pinery_temp_annual, aes(x = year, y = mean_temp), colour="lightsalmon4") +
#Aesthetics
xlab("Year") +
ylab("Mean Annual Temperature (°C)") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
#ggtitle("Average Annual Temperature in \n Southwestern Ontario from 1883 to 2017") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
annual_temp_plot## Warning: Removed 1 row(s) containing missing values (geom_path).
#ggsave('annual_temp.pdf', annual_temp_plot, units = 'cm', width = 15, height = 10)
annual_precip_plot <- ggplot() +
# London Data
geom_line(data=london_precip_annual, aes(x = year, y = total_precip), colour = "lightseagreen", linetype = "dashed", alpha = 0.8) +
# geom_smooth(data=london_precip_annual, aes(x = year, y = total_precip), se = TRUE,
# method = 'lm', size = 0.8, colour = "black", fill = "blue") +
# Pinery Data
geom_line(data=pinery_precip_annual, aes(x = year, y = total_precip, group = 1), colour="lightsalmon4")+
# geom_smooth(data=pinery_precip_annual, aes(x = year, y = total_precip), se = TRUE,
# method = 'lm', size = 0.8, colour = "black", fill = "darkgreen") +
#Aesthetics
xlab("Year") +
ylab("Total Annual Precipitation (mm)") +
scale_x_continuous(breaks=seq(1880,2020,20)) +
#ggtitle("Annual Precipitation in Southwestern Ontario \nfrom 1883 to 2017") +
scale_colour_manual(name="Location", breaks = c("London", "Pinery"),
values = c("blue","darkgreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
annual_precip_plot## Warning: Removed 17 row(s) containing missing values (geom_path).
#ggsave('annual_precip.pdf', annual_precip_plot, units = 'cm', width = 15, height = 10)
annual_plot <- grid.arrange(annual_temp_plot, annual_precip_plot, ncol= 2)## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 17 row(s) containing missing values (geom_path).
#ggsave('annual_plot.pdf', annual_plot, units = 'cm', width = 23, height = 8)6 Seasonal Data
winter_t <- pinery_temp %>% filter(month == 1 | month == 2 | month == 3)
winter_t$season <- "winter"
spring_t <- pinery_temp %>% filter(month == 4 | month == 5 | month == 6)
spring_t$season <- "spring"
summer_t <- pinery_temp %>% filter(month == 7 | month == 8 | month == 9)
summer_t$season <- "summer"
fall_t <- pinery_temp %>% filter(month == 10 | month == 11 | month == 12)
fall_t$season <- "fall"
pinery_temp <- rbind(winter_t, spring_t, summer_t, fall_t)
pinery_temp$season <- as.factor(pinery_temp$season) #Make a factor
str(pinery_temp)## 'data.frame': 1944 obs. of 5 variables:
## $ year : int 1885 1885 1885 1886 1886 1886 1887 1887 1887 1888 ...
## $ month : int 1 2 3 1 2 3 1 2 3 1 ...
## $ mean_temp: num -9.4 -13.5 -8.9 -7.2 -6.9 -1 -8.4 -5.4 -3.8 -9.6 ...
## $ station : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ season : Factor w/ 4 levels "fall","spring",..: 4 4 4 4 4 4 4 4 4 4 ...pinery_temp_station <- pinery_temp %>%
group_by(year, season, station) %>%
summarise(mean_temp = mean(mean_temp))## `summarise()` regrouping output by 'year', 'season' (override with `.groups` argument)pinery_temp_season <- pinery_temp_station %>%
group_by(year, season) %>%
summarise(mean_temp_season = mean(mean_temp) %>% round(2), sd_temp = sd(mean_temp) %>% round(2), n = n()) %>%
mutate(se_temp = (sd_temp / sqrt(n)) %>% round(2), lower = mean_temp_season - se_temp, upper = mean_temp_season + se_temp) ## `summarise()` regrouping output by 'year' (override with `.groups` argument)winter_p <- pinery_precip %>% filter(month == 1 | month == 2 | month == 3)
winter_p$season <- "winter"
spring_p <- pinery_precip %>% filter(month == 4 | month == 5 | month == 6)
spring_p$season <- "spring"
summer_p <- pinery_precip %>% filter(month == 7 | month == 8 | month == 9)
summer_p$season <- "summer"
fall_p <- pinery_precip %>% filter(month == 10 | month == 11 | month == 12)
fall_p$season <- "fall"
pinery_precip <- rbind(winter_p, spring_p, summer_p, fall_p)
pinery_precip$season <- as.factor(pinery_precip$season) #Make a factor
pinery_precip_station <- pinery_precip %>%
group_by(year, season, station) %>%
summarise(total_precip = mean(total_precip))## `summarise()` regrouping output by 'year', 'season' (override with `.groups` argument)pinery_precip_season <- pinery_precip_station %>%
group_by(year, season) %>%
summarise(total_precip_season = mean(total_precip) %>% round(2), sd_precip = sd(total_precip) %>% round(2), n = n()) %>%
mutate(se_precip = sd_precip / sqrt(n), lower = total_precip_season - sd_precip, upper = total_precip_season + sd_precip %>% round(2)) ## `summarise()` regrouping output by 'year' (override with `.groups` argument)winter_tl <- london_temp %>% filter(month == "Jan" | month == "Feb" | month == "Mar")
winter_tl$season <- "winter"
spring_tl <- london_temp %>% filter(month == "Apr" | month == "May" | month == "Jun")
spring_tl$season <- "spring"
summer_tl <- london_temp %>% filter(month == "Jul" | month == "Aug" | month == "Sep")
summer_tl$season <- "summer"
fall_tl <- london_temp %>% filter(month == "Oct" | month == "Nov" | month == "Dec")
fall_tl$season <- "fall"
london_temp <- rbind(winter_tl, spring_tl, summer_tl, fall_tl)
london_temp$season <- as.factor(london_temp$season) #Make a factor
london_temp_season <- london_temp %>%
group_by(year, season) %>%
summarise(mean_temp_season = mean(mean_temp))## `summarise()` regrouping output by 'year' (override with `.groups` argument)winter_pl <- london_precip %>% filter(month == "Jan" | month == "Feb" | month == "Mar")
winter_pl$season <- "winter"
spring_pl <- london_precip %>% filter(month == "Apr" | month == "May" | month == "Jun")
spring_pl$season <- "spring"
summer_pl <- london_precip %>% filter(month == "Jul" | month == "Aug" | month == "Sep")
summer_pl$season <- "summer"
fall_pl <- london_precip %>% filter(month == "Oct" | month == "Nov" | month == "Dec")
fall_pl$season <- "fall"
london_precip <- rbind(winter_pl, spring_pl, summer_pl, fall_pl)
london_precip$season <- as.factor(london_precip$season) #Make a factor
london_precip_season <- london_precip %>%
group_by(year, season) %>%
summarise(total_precip_season = mean(total_precip))## `summarise()` regrouping output by 'year' (override with `.groups` argument)london_temp_season$season <- factor(london_temp_season$season, levels = c("winter", "spring", "summer", "fall"))
pinery_temp_season$season <- factor(pinery_temp_season$season, levels = c("winter", "spring", "summer", "fall"))pinery_temp_season$location <- "Pinery"
london_temp_season$location <- "London"
season_temp <- rbind(pinery_temp_season, london_temp_season)
pinery_precip_season$location <- "Pinery"
london_precip_season$location <- "London"
season_precip <- rbind(pinery_precip_season, london_precip_season)
season_temp$location <- factor(season_temp$location, levels = c("Pinery", "London"))
season_precip$location <- factor(season_precip$location, levels = c("Pinery", "London"))6.1 Temperature
w_temp <- ggplot() +
geom_line(data=(season_temp %>% filter(season == "winter")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Temperature (°C)") + ggtitle("Winter") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
sp_temp <- ggplot() +
geom_line(data=(season_temp %>% filter(season == "spring")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Temperature (°C)") + ggtitle("Spring") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
su_temp <- ggplot() +
geom_line(data=(season_temp %>% filter(season == "summer")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Temperature (°C)") + ggtitle("Summer") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
fa_temp <- ggplot() +
geom_line(data=(season_temp %>% filter(season == "fall")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Temperature (°C)") + ggtitle("Fall") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
seasonal_temp_plot <- grid.arrange(w_temp, sp_temp, su_temp, fa_temp, nrow = 2, ncol = 2)## Warning: Removed 1 row(s) containing missing values (geom_path).
#ggsave('seasonal_temperature.pdf', seasonal_temp_plot, units = 'cm', width = 20, height = 12)6.2 Precipitation
w_precip <- ggplot() +
geom_line(data=(season_precip %>% filter(season == "winter")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Winter") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
sp_precip <- ggplot() +
geom_line(data=(season_precip %>% filter(season == "spring")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Spring") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
su_precip <- ggplot() +
geom_line(data=(season_precip %>% filter(season == "summer")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Summer") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
fa_precip <- ggplot() +
geom_line(data=(season_precip %>% filter(season == "fall")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Fall") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
seasonal_precip_plot <- grid.arrange(w_precip, sp_precip, su_precip, fa_precip, nrow = 2, ncol = 2)## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 15 row(s) containing missing values (geom_path).
## Warning: Removed 15 row(s) containing missing values (geom_path).## Warning: Removed 16 row(s) containing missing values (geom_path).
#ggsave('seasonal_precipitation.pdf', seasonal_precip_plot, units = 'cm', width = 20, height = 12)6.3 Temperature Anomaly
https://tamino.wordpress.com/2018/07/29/why-use-temperature-anomaly/
Average defined as 1951-1980 (NASA standard)
pinery_average <- season_temp %>% filter(location == "Pinery") %>% filter(year > 1950 & year < 1981)
pinery_average <- pinery_average %>% group_by(season) %>% summarise(pinery_average = mean(mean_temp_season))## `summarise()` ungrouping output (override with `.groups` argument)london_average <- season_temp %>% filter(location == "London") %>% filter(year > 1950 & year < 1981)
london_average <- london_average %>% group_by(season) %>% summarise(london_average = mean(mean_temp_season))## `summarise()` ungrouping output (override with `.groups` argument)pinery_anomaly <- right_join(season_temp, pinery_average)## Joining, by = "season"pinery_anomaly <- mutate(pinery_anomaly, anomaly = mean_temp_season - pinery_average) %>% filter(location == "Pinery")
london_anomaly <- right_join(season_temp, london_average)## Joining, by = "season"london_anomaly <- mutate(london_anomaly, anomaly = mean_temp_season - london_average) %>% filter(location == "London")
season_anomaly <- full_join(pinery_anomaly, london_anomaly) %>% select(-pinery_average, -london_average)## Joining, by = c("year", "season", "mean_temp_season", "sd_temp", "n", "se_temp", "lower", "upper", "location", "anomaly")w_anomaly <- ggplot() +
geom_line(data=(season_anomaly %>% filter(season == "winter")), aes(x = year, y = anomaly, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Winter") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
sp_anomaly <- ggplot() +
geom_line(data=(season_anomaly %>% filter(season == "spring")), aes(x = year, y = anomaly, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Spring") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
su_anomaly <- ggplot() +
geom_line(data=(season_anomaly %>% filter(season == "summer")), aes(x = year, y = anomaly, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Summer") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
fa_anomaly <- ggplot() +
geom_line(data=(season_anomaly %>% filter(season == "fall")), aes(x = year, y = anomaly, col = location, linetype = location)) +
#Aesthetics
xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Fall") +
scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
seasonal_anomaly_plot <- grid.arrange(w_anomaly, sp_anomaly, su_anomaly, fa_anomaly, nrow = 2, ncol = 2)## Warning: Removed 1 row(s) containing missing values (geom_path).
#ggsave('seasonal_anomaly.pdf', seasonal_anomaly_plot, units = 'cm', width = 20, height = 12)7 Save Outputs
write.csv(PL.clean, "outputs/full_climate.csv")
write.csv(PL.temp4, "outputs/full_temp.csv")
write.csv(PL.precip4, "outputs/full_precip.csv")
write.csv(pinery_temp_mean, "outputs/pinery_monthly_temperature.csv")
write.csv(pinery_precip_mean, "outputs/pinery_monthly_precipitation.csv")
write.csv(season_temp, file = "outputs/pinery_season_temperature.csv")
write.csv(season_precip, file = "outputs/pinery_season_precipitation.csv")
write.csv(pinery_temp_annual, file = "outputs/pinery_annual_temp.csv")8 Reproducibility
Sys.time()## [1] "2020-08-06 16:25:48 PDT"git2r::repository()## Local: master /Users/JenBaron/Documents/UWO/UWO NSERC/Statistical Analysis/Climate/Climate_Pinery
## Remote: master @ origin (https://github.com/JenBaron/Climate_Pinery.git)
## Head: [78f8528] 2020-08-05: filled .rmdsessionInfo()## R version 4.0.0 (2020-04-24)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gridExtra_2.3 GGally_1.5.0 forcats_0.5.0 stringr_1.4.0
## [5] dplyr_1.0.0 purrr_0.3.4 readr_1.3.1 tidyr_1.0.2
## [9] tibble_3.0.1 tidyverse_1.3.0 ggplot2_3.3.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 xfun_0.14 haven_2.2.0 lattice_0.20-41
## [5] colorspace_1.4-1 vctrs_0.3.1 generics_0.0.2 htmltools_0.4.0
## [9] yaml_2.2.1 utf8_1.1.4 rlang_0.4.6 pillar_1.4.4
## [13] glue_1.4.1 withr_2.2.0 DBI_1.1.0 RColorBrewer_1.1-2
## [17] dbplyr_1.4.3 modelr_0.1.6 readxl_1.3.1 plyr_1.8.6
## [21] lifecycle_0.2.0 munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0
## [25] rvest_0.3.5 evaluate_0.14 labeling_0.3 knitr_1.28
## [29] fansi_0.4.1 broom_0.5.6 Rcpp_1.0.4.6 scales_1.1.1
## [33] backports_1.1.7 jsonlite_1.6.1 farver_2.0.3 fs_1.4.1
## [37] hms_0.5.3 digest_0.6.25 stringi_1.4.6 grid_4.0.0
## [41] cli_2.0.2 tools_4.0.0 magrittr_1.5 crayon_1.3.4
## [45] pkgconfig_2.0.3 ellipsis_0.3.1 xml2_1.3.2 reprex_0.3.0
## [49] lubridate_1.7.8 reshape_0.8.8 assertthat_0.2.1 rmarkdown_2.1
## [53] httr_1.4.1 rstudioapi_0.11 R6_2.4.1 git2r_0.26.1
## [57] nlme_3.1-147 compiler_4.0.0